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Abstract 



It has recently been shown that the Helmholtz free energy difference be- 
tween two equihbrium configurations of a system may be obtained from an 
ensemble oi finite-time (nonequilibrium) measurements of the work performed 
in switching an external parameter of the system. Here this result is estab- 
lished, as an identity, within the master equation formalism. Examples are 
discussed and numerical illustrations provided. 
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INTRODUCTION 



Consider some finite classical system which depends on an external parameter, A. For 
instance, the system might be a lattice of coupled classical spins, and A may denote the 
strength of an externally applied magnetic field; or, the system might be a gas of particles, 
and A a parameter specifying the volume of a box enclosing the gas. Now suppose that, 
after allowing the system to come to equilibrium with a heat reservoir at temperature T, we 
change, or "switch" , the external parameter, infinitely slowly, from an initial (say, A = 0) to 
a final (A = 1) value. The system will remain in quasi-static equilibrium with the reservoir 
throughout the switching process, and the total work performed on the system will equal 
the Helmholtz free energy difference between the initial and final configurations |l|]: 

W^ = AF = Fi- Fo. (1) 

Here Fx denotes the equilibrium free energy of the system at temperature T, for a fixed value 
of A. The subscript on W reminds us that this result is valid for infinitely slow switching of 
the parameter. 

Now, what happens if, after allowing the system and reservoir to equilibrate, we switch 
the value of A at a finite rate? In this case the system will lag behind quasi-static equilibrium 
with the reservoir, and the total work will depend on the microscopic initial conditions of 
system and reservoir. Thus, an ensemble of such switching measurements — each prepared 
by first allowing the system to equilibrate with the reservoir — will yield a distribution of 
values of W. Let p{W,ts) denote this distribution, where tg is the "switching time" over 
which the value of A is changed from to 1. (Without loss of generality, assume a uniform 
switching rate: A = t~^-) In other words, p{W,ts) dW is the probability that the work 
performed in switching A from to 1, over a time tg, will fall between W and W + dW . In 
the limit ts —>■ oo, we get W = AF (Eq.|I|), and so p —>■ 6{W — AF) in this limit. For finite 
tg, however: (1) the distribution p acquires a finite width, reflecting the fluctuations in W 
from one switching measurement to the next; and (2) its centroid shifts to the right: 

W = J dW p{W,ts)W > AF, (2) 

result of dissipation (see Fig.|l]) [@,^. 
Eq.|l] gives the free energy difference AF in terms of the work performed during a single 
infinite-time (quasi-static) process. By contrast, Eq.|^ gives an upper bound on AF, from 
an ensemble of finite-time (hence, nonequilibrium) repetitions of the switching process. Re- 
cently 0, it was shown that one can in fact extract the value of AF itself, not just an 
upper bound, from the information contained in p{W,ts). Specifically, the following result 
was shown to be valid for any switching time tg-. 

exp -(3W = JdW p{W, ts) exp -(3W = exp -(3AF, (3) 

where (3~^ = ksT. This result gives the value of an equilibrium quantity, AF, in terms of 
an ensemble of finite-time, nonequilibrium measurements. As discussed in Ref. and in 
Section ^ below, the inequality W > AF, Eq.^ follows immediately from Eq.^. 
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Eq.^ was derived by treating the system of interest and reservoir, coupled together, as 
a large, isolated, classical Hamiltonian system. The Hamiltonian governing the motion in 
the full phase space was taken to be a sum of three terms: one for the system of interest 
(Hx), one for the reservoir, and a final term, hint, coupling the two. The magnitude of 
the interaction term hint was explicitly assumed to be negligible in comparison with the 
other two terms. With these assumptions, Eq.^ follows from the properties of Hamiltonian 
evolution (see Ref. [Q for details). 

The purpose of the present paper is to re-derive the same result within a different frame- 
work. Instead of considering deterministic evolution in the full phase space — containing 
both the system of interest and reservoir degrees of freedom — we will treat only the evo- 
lution of the system of interest itself, described by a trajectory z(t). This evolution will be 
stochastic rather than deterministic, and will be governed by a master equation. We will 
assume that this stochastic evolution is Markovian, and that it satisfies detailed balance. 

In presenting this alternative derivation, we are motivated by several factors. First, 
master equations are a common tool in statistical physics, therefore it is reassuring to see 
that Eq.^ follows in a natural way from the master equation approach. Secondly, in this 
derivation there is no need to explicitly assume weak coupling between system and reservoir, 
since a reservoir per se does not enter into the analysis; given the assumptions stated below, 
Eq.|^ is identically true. Finally, one might come away from Ref. Q with the feeling that that 
the validity of Eq.|^ depends directly on the properties of Hamiltonian evolution, in particular 
Liouville's theorem. The treatment herein dispels this notion: Hamilton's equations appear 
nowhere in the derivation. This point is particularly relevant in the context of numerical 
simulations, where the evolution of a thermostatted system is often realized with the use of 
non-Hamiltonian equations of motion. 

The plan of this paper is as follows. In Section |I| we establish notation and terminology, 
and specify precisely our assumptions regarding Markovian evolution and detailed balance. 
In Section |T| we derive our central result. In Section JTT| we consider specific examples 
of stochastic processes for which the result is valid. In Section |I^ we briefly discuss the 
possible utility of Eq.§ to the numerical computation of free energy differences. In Section 
|V] we illustrate our central result with numerical simulations, and we conclude with a few 
remarks in Section W^. Two appendices provide derivations of results used in the main text. 



I. PRELIMINARIES 

We begin by specifying precisely what we mean by the terms work and free energy 
throughout this paper. We assume that there exists a phase space of variables (e.g. the 
positions and momenta of constituent particles), such that the instantaneous microscopic 
state of our system is completely described by the values of these variables. Let z denote 
a point in this phase space. The evolution of our system with time is then described by a 
trajectory z{t). The kind of evolution which will most interest us is not the evolution of 
an isolated system, but rather that of a system in contact with a heat bath. Hence the 
trajectory z{t) will in general be stochastic, reflecting the "random" influence of the heat 
bath. 

Next, assume the existence of a parameter-dependent Hamiltonian Hx{z). The Hamil- 
tonian is just a function which, for a fixed value of A, gives the total energy of the system 
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of interest, in terms of its instantaneous state (z). The value of A thus parametrizes the 
external forces acting on the system (arising from, e.g., external fields, confining potentials, 
etc.) When the system is isolated, Hx happens to generate the time evolution of the system, 
through Hamilton's equations, but we will not make use of this property in deriving our 
central result. 

Given Hx{z), we now define 

Zx{f3)= Jdzexp-f3Hx{z) (4a) 
F,(/?)^-/?-MnZ,(/?), (4b) 

where /3 is a real, positive constant. (We are being a bit cavalier with units here: Z\ should 
be divided by a constant to make it dimensionless. However, since this constant only shifts 
the value of Fx by a fixed amount, and therefore does not affect the free energy difference 
Fi — Fq, we ignore it.) Z\{I3) and Fx{j3) are of course the partition function and free energy 
of the system, respectively, and is the temperature, in units of energy (/5 = l/ksT). 
However, in following the derivation presented below, it may be most convenient to view 
these quantities abstractly, rather than in connection with their physical significance: f3 
is just some positive constant, and Z\{[3) and F\{[3) are the functions of j3 defined by 
Eq.^. In what follows, we will never compare free energies or partition functions at different 
temperatures, hence the dependence of Fx and Zx on {3 will be left implicit. 

The central quantity of interest in this paper will be the free energy difference AF = 
Fi — -Fq, for a fixed value of (3. 

As mentioned, the time evolution of the system of interest is described entirely by a 
(stochastic) phase space trajectory z(t). In general, this evolution will depend on the ex- 
ternally imposed time- dependence of A, describing the changing external fields to which the 
system is subject. Let us consider the evolution of the system from an initial time, t = 0, 
to a final time, t = is, over which the value of A is "switched" from to 1 at a uniform 
rate: A(t) = t/ts- Given this time- dependence of A, and given the trajectory z(t) describing 
the evolution of the system, the total work performed on the system is the time integral of 
XdHx/dX along the trajectory: 

where X = dX/dt = t~^. For the evolution of an isolated Hamiltonian system, this reduces to 
W = Hi{7j{ts)) — Ho{z{0)), by Hamilton's equations; in this case the work performed on the 
system is just the change in its energy. For a system in contact with a heat bath, however, 
this no longer holds, since there is a constant exchange of energy with the bath. 

We will assume that the evolution of our system in phase space is a Markov process 
This means that the stochastic evolution z{t) is entirely characterized by the transition 
probability function P{z', t\z, t + At). This gives the probability distribution for finding the 
system in a state z at time t + At, given that at an earlier time t it was known to be at z' . 
Taking the derivative of P with respect to At and evaluating at At — > O"*", we get a function 

(6) 



i?(z',z;t) 



d 

d{At) 



P(z',t|z,t + At) 
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which gives the instantaneous transition rate from z' to z, at time t. The dependence of R on 
time arises through whatever external parameters of the system and reservoir are available 
and time-dependent. In our case we assume only one such parameter, A — characterizing 
external forces — therefore we write 



R{z',z;t)^Rxiz',z). (7) 

In other words, the instantaneous transition rate R from z' to z depends on t only through 
the value of A(t). 

Let us now shift our focus from the description of a single system, to that of an ensemble 
of systems, each evolving according to the stochastic Markov process just described. This 
ensemble represents infinitely many independent realizations of the switching process. If 
f{z,t) denotes the time-dependent distribution of this ensemble in phase space, then this 
distribution obeys 



9f^ 



z,t)= Jdz'f{z',t)R^{z',z), (8a) 



dt 

where A = A(t). We will abbreviate this as 

df 

^ = /, (8b) 

where Rx is a linear operator acting on the space of phase space densities /. Eq.|^ is our 
master equation; when the time-dependence of A is known, and an initial distribution /o is 
specified, Eq.H uniquely determines the subsequent evolution of the phase space density /. 

In addition to the Markov assumption (Eq.^, we will impose another assumption on our 
stochastic process: detailed balance. If A is held fixed, z(t) becomes a stationary Markov 
process, describing the evolution of a system in contact with a heat reservoir, when the exter- 
nal forces acting on the system are time-independent. Under such evolution, the canonical 
distribution in phase space (corresponding to the fixed value of A) ought to be invariant. By 
Eq.|^, this is equivalent to the statement that Rx annihilates the canonical distribution: 

Rxexp-pHxiz) = 0. (9) 

This places a condition on the linear operator Rx- We assume that our stochastic process 
satisfies this condition, and will refer to this assumption as detailed balance. (Usually, 
detailed balance is expressed as the following, somewhat stronger assumuption: 

R\{z',z) ^ exp -(3Hx{z) 

Rx{z,z') exp-pHx{z')- ^ ^ 

Eq.p follows immediately from Eq.|lO|: just set the products of the cross terms equal, then 
integrate over z'. The distinction between Eqs.^ and |10| is of little importance in the present 
context, so for simplicity we refer to Eq.^ as detailed balance.) 

Having been led to assume detailed balance by considering the behavior of the system 
when A is fixed, it may seem natural to make another assumption. Namely, if our stochastic 
process is meant to describe a system in contact with a reservoir, then for fixed A we expect 
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the system to thermalize: after an initial relaxation time, the system samples its phase 
space canonically; equivalently, an arbitrary initial ensemble /o(z) relaxes to a canonical 
distribution in phase space. We may express this assumption formally, as: 

hm f/A(t)/o(z) = ^exp-/?i7,(z), (11) 

for any normalized /o(z). The operator tf\{t) = expRxt appearing on the left is just the 
evolution operator corresponding to the equation of motion df/dt = R\f, for fixed A. 
Eq.|TT| says that any initial distribution /o(z) relaxes to a canonical distribution, and then 
stays there. We will refer to this as the assumption of thermalization. Note that while 
thermalization (Eq.|lT]) implies detailed balance (Eq.|D, the converse is not true. Since it 
will turn out that the proof of EqJ^ requires only the weaker assumption of detailed balance, 
and not the stronger thermalization assumption, we will assume in what follows that Eq.|^ 
holds, but not necessarily Eq.^^ unless explicitly stated. 

[Physically, we expect both thermalization and detailed balance to hold, as long as our 
master equation describes a system in contact with a genuine heat reservoir. However, one 
can easily imagine situations in which the evolution satisfies Eq.^ but not Eq.jT^. A specific 
example, which we will analyze in Section |Tl|, is that of an isolated Hamiltonian system. 
Note that if thermalization does not hold, then the validity of Eq.|I|, Woo = ^F, is not 
guaranteed, since that result assumes that during an infinitely slow switching process the 
system remains in quasi-static equilibrium with a reservoir. Eq.|] nevertheless remains valid, 
provided detailed balance is satisfied.] 

At this point we are ready to proceed with a proof of EqJ^. Concrete examples of 
stochastic processes satisfying the assumptions made in this section will be discussed in 
Section |T|, and numerically simulated in Section M. 



II. DERIVATION 

The overbar appearing on the left side of Eq.|^ denotes an average over an infinite en- 
semble of independent realizations (repetitions) of the switching process. Each realization 
is described by a trajectory z(t), < t < tg, specifying the evolution of the complete set 
of phase space variables as the external parameter A is switched from to 1. The entire 
ensemble is then described by a time-dependent phase space density: at any time t, f(z,t) 
represents a snapshot of the distribution of trajectories in phase space. Since we assume 
that the system equilibrates with the reservoir (with A held fixed at 0) before the start of 
each realization, we have a canonical distribution of initial conditions 0: 

f{z,0) = Z^' exp~PHo{z). (12) 

During the switching process, however, the ensemble does not (in general) remain in instan- 
taneous canonical equilibrium. In other words, although the distribution / = Z^^ exp —13H\ 
is a (stationary) solution of Eq.^ for A fixed, it is not, in general, a solution of Eq.^ when A 
depends on time. Thus, for t > 0, a snapshot of the trajectories will reveal a distribution 
/(z,t) which lags behind the canonical distribution corresponding to \{t). The amount of 
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lag present by the time a given value of A is reached, will depend on how rapidly or slowly 
we performed the switching on the way to that value. 

For every trajectory z{t) in our ensemble, we can compute the total work W performed 
on the system (Eq.|^). Our task is now to evaluate the ensemble average of exp —(3W. To 
do this, let us first define, for a given trajectory z(t), a function w{t) which is the "work 
accumulated" up to time t: 

w{t) = l^dt'X^{z{t')). (13) 

Thus, W = w{ts)- Now consider all those trajectories in the ensemble which happen to 
pass through the phase space point z at time t, and let Q{z,t) denote the average value of 
exp —l3w{t), over this particular subset of trajectories. Finally, define 

g{z,t) = f{z,t)Qiz,t). (14) 

Note that g{z, 0) = /(z, 0), since w{0) = for every trajectory. Given these definitions, the 
ensemble average of exp —(3W may be expressed as: 

exp-/?Vr = j dzg{z,ts). (15) 

We will now solve for g{z,ts). 
The function g{z,t) obeys 

with A = A(t). To see this, imagine for a moment that each trajectory z{t) in the ensemble 
represents a "particle" moving about in phase space. Furthermore, assume that the "mass" 
of each particle is time-dependent, and is given by fi{t) = exp —(3w{t). Each particle thus 
begins with mass unity: /i(0) = 1. The function Q{z,t) is then the average mass of those 
particles which are found at a point z at time t, and g{z,t) represents the time-dependent 
mass density in phase space (after normalization of the number density to unity: J f dz = 1). 
This mass density is time-dependent for two reasons. First, (A) the mass of each particle 
changes with time, at a rate proportional to its instantaneous mass: 

m = -Pw{t)fi{t) = -/5A^(z(t))/x(t). (17) 
This contributes a term 

: -pX^{z) g{z,t) (18) 

to dg/dt. Secondly, (B) the mass density evolves due to the flow of particles, described by 
our master equation, df jdt = R\ f. This flow of particles contributes a term 

=fl.9. (19) 
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Adding these contributions, A and B, gives Eq.|T6[ An alternative derivation of this evolution 
equation, based on a path integral formulation, is given in Appendix A. 

Given the initial conditions (^(z, 0) = /(z, 0) = Zq^ exp —PHq{z), Eq.|16| is solved by 

giz,t) = Z^'exp-pH.iz) = |^/f (z), (20) 

where A = A(t), and denotes the canonical distribution in phase space, for a given value 
of A. This result is easily verified with the help of Eq.|^. Then, using Eq.[l^, we finally arrive 

at 

exp -PW = Zq^ [ dz exp -(3Hi{z) = |i = exp -/3AF. (21) 

Q.E.D. 

A different proof of Eq.^ has been discovered by Gavin E. Crooks ||^ . 

It is worthwhile to draw attention to a curious feature of the evolution of our "mass 
density", g{z,t). Recall that the evolution of f{z,t) depends nontrivially on the rate at 
which we switch A: at finite switching rates, /(z, t) lags behind the instantaneous equilibrium 
distribution. By contrast, the time dependence of g{z,t) is very simple: the mass density 
is determined uniquely by the instantaneous value of A; see Eq.^ Thus, no matter how 
slowly or rapidly we switch A from to 1, g{z, t) will always evolve through exactly the same 
continuous sequence of "canonical" mass densities, specified by Eq.|20|; g never develops a 
lag, in the sense that / does. 

Eq.^ also implies the result 

exp -/3wit) = 1^ = exp - Fq), (22) 

A = A(t), which is in effect a restatement of our central result, EqJ^. 

III. EXAMPLES 

We have shown that Eq.^ is true for Markovian processes satisfying detailed balance. 
Let us consider a few examples of such processes. 

A. Hamiltonian evolution 

As a first case, we take deterministic evolution under the Hamiltonian Hx, as A is varied 
from to 1: 

z = {z,H,}, (23) 

where {■,■} denotes the Poisson bracket. Ordinarily, one would not call this a stochastic 
process, but of course it may be viewed as a special case of such. This evolution is Markovian, 
and Rx is just the Poisson bracket operator: 

Rxf = {HxJ}. (24) 
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It immediately follows that detailed balance (Eq.^ is satisfied. The central result of this 
paper then tells us that, if we (1) start with a canonical distribution of initial conditions, 
/(z, 0) = Zq^ exp —PHq{z), and (2) allow trajectories to evolve deterministically from these 
initial conditions (Eq. |23|) as A is varied from to 1, then this ensemble of trajectories will 
satisfy Eq.^, regardless of how slowly or quickly we perform the switching. This result was 
proven more directly in Ref. 0]. 

Since evolution under Eq.^ describes an isolated system (no heat bath), the thermaliza- 
tion assumption (Eq.[Tl|) is not met, and the work performed in the limit of infinitely slow 
switching will in general not equal the free energy difference: W^o 7^ AF. Eq.|^ nevertheless 
remains true, even in this limit. 

This last statement is easily illustrated with a one- dimensional harmonic oscillator. Let 
the Hamiltonian be 

Hx{x,p) = ^ + u;l^. (25) 

For a given temperature, the partition function is given by Zx = 27c/i3uJx, and the free energy 
difference is 

AF = -r'ln|^ = r'ln-. (26) 

Zq ujo 

Let us now imagine a trajectory z{t) evolving under this Hamiltonian, as A is changed 
infinitely slowly from to 1; assume for specificity that ui > cjq- Since H\/uj\ is an 
adiabatic invariant for the harmonic oscillator p|, we have Ei = {uji/uq) Eq, where Eq [Ei) 
is the initial (final) energy of the oscillator. The work performed on an isolated system is 
equal to the change in its energy, so we get Woo = [(t^i/t^o) — l]-E'o- A canonical distribution 
of initial energies Eq then leads, after some algebra, to the following ensemble distribution 
of values of W: 

lim piW,Q = ^^expf-^^)^(l^), (27) 

where 6 denotes the unit step function. It is straighforward to verify that this distribution 
satisfies Eq.|^. 

As another example, consider a single particle bouncing around inside a three- 
dimensional cavity with hard walls, where the shape of the cavity is a function of A. Let 
Va denote the volume of the cavity. The free energy difference is then AF = f3~^ ln(Vo/Vi). 
Assume that, when A is held fixed, the motion of the particle is ergodic over the energy shell 
(the 5-dimensional surface of constant energy in 6-dimensional phase space); also assume 
that Vi < Vq. Now, imagine that we allow the particle to evolve as A is switched infinitely 
slowly from to 1. For this system the quantity H^^Vx is an adiabatic invariant 0, therefore 
the total work performed on a particle with initial energy Eq is: Woo = [(Vo/Vi)^'''^ — l]Eo. 
Taking an ensemble of such particles, defined by a canonical distribution of initial conditions 
(at A = 0), we get 

hm p{W,Q = ( ' expf-^V(^)' (28) 



where r = (Vq/Vi)^/^ — 1. As with the example of the harmonic oscillator, it is straightfor- 
ward to show that this distribution p satisfies Eq.^. 
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B. Langevin evolution 



Next, let us consider modifying Hamilton's equations, by adding both a frictional and a 
stochastic force: 

x = p (29a) 
p=-d,Vx-w + F{t). (29b) 

We have assumed a one-degree-of-freedom system, and a Hamiltonian Hx = + V\{x). 
Let us furthermore take the stochastic force, F{t), to represent white noise, with an auto- 
correlation function given by 

{F{t,)F{t,)) = Dp6{h-t,). (30) 

Finally, let us impose a fluctuation-dissipation relation between the frictional and stochastic 
forces: 

7 = PDp/2. (31) 

An ensemble of trajectories evolving under the stochastic process just defined satisfies the 
Fokker-Planck equation, 

+ (32, 

Since Eq.^ is of the form df /dt = Rxf, this process is Markovian. Furthermore, inspection 
reveals that Rx satisfies Eq.^, i.e. detailed balance holds. Thus, the conditions for the 
validity of Eq.§ are met: given an ensemble of systems evolving under this stochastic process, 
as A is changed from to 1 over a time tg, and given an initial distribution /(z,0) = 
Zq^ exp —PHo(z), we are guaranteed that this ensemble will satisfy exp —(3W = exp — 
for any switching time tg- 

Eq.^is a Langevin set of equations. If A is held fixed, an ensemble of trajectories evolving 
under these stochastic equations of motion will relax to a canonical distribution on phase 



space. In other words, the evolution specified by Eqs.|^ to ^ satisfies not only detailed 
balance, but the stronger assumption of thermalization as well. If we start with a canonical 
distribution, at A = 0, and then switch the value of the external parameter infinitely slowly 
to A = 1, then the ensemble of trajectories will remain in quasi-static equilibrium over the 
course of the switching process: 

f{z,t) = Z^'exp-l3Hxiz) , X = X{t) , ^ oo. (33) 

The creeping value of A thus drags f{z,t) through a continuous sequence of canonical dis- 
tributions. Furthermore, in this limit, the "work accumulated" will be the same function 
of time for every trajectory in the ensemble: w{t) = Fx — Fq; each trajectory represents a 
system evolving in quasi-static equilibrium with the reservoir. Thus, the quantity Q{z,t) 
defined earlier is given, in this quasi-static limit, by 

Q(z, t) = exp -(3{Fx - Fq) = Zx/Z^. (34) 
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Zq ^ exp —PH, 



(Eq.| 



It then follows that the function g = fQ (Eq.|^) satisfies g{z, t) 
from which we get exp —f3W = exp —f3AF. Of course, when we switch A from to 1 over a 
finite time tg, neither Eq.|33| nor Eq.|33 will in general hold; nevertheless, Eq.EOl and therefore 



A 



our central result, will remain valid, since Eq.|16| governs the time evolution of g regardless 
of how slowly or quickly we switch the external parameter. 



C. Isothermal molecular dynamics 

The Langevin equations above provide a simple method for numerically simulating the 
evolution of a thermostatted system (i.e. a system in contact with a heat bath), without 
explicitly simulating the many degrees of freedom of the bath. The term F{t) may be imple- 
mented by generating a small, random momentum "kick" at each time step in the numerical 
integration of the equations of motion. This method works both for static Hamiltonians and 
— as in the case considered in the present paper — for Hamiltonians made time- dependent 
through the variation of an external parameter. In the past decade or so, methods for 
using explicitly deterministic equations of motion to simulate thermostatted systems, have 
proven very useful |T^. These generally go under the name of isothermal molecular dy- 



namics (IMD). Typically, the heat bath is represented by one or more "extra" degrees of 
freedom. Then, in the extended phase space which includes both the system of interest 
and the extra degree (s) of freedom, the evolution is governed by a set of deterministic (but 
non-Hamiltonian) equations of motion. These are tailored so that, when the Hamiltonian 
describing the system of interest is static, the variables representing the system of interest 
explore phase space canonically, at least to a good approximation. 

An example of an IMD scheme is Nose-Hoover dynamics |TT||T^, represented by the 
following equations of motion: 

{<i = P , P = -VVa - Cp}^^ (35a) 
C = (K/Ko - \)It\ (35b) 

We have again assumed a kinetic-plus-potential Hamiltonian; the index n runs over all D 
degrees of freedom of the system, K = J^nPnf^ is the total kinetic energy of the system, 
Kq = D/2f3 is the thermal average of K, and the parameter r acts as a relaxation time. 
Let us imagine a trajectory evolving in the extended phase space, (z, ^)-space, under these 
equations of motion, as A is switched from to 1. The work W performed on the system 
of interest is defined, as before, by Eq.^. While this expression does not explicitly contain 
the bath variable C, W is nevertheless a function of the full set of initial conditions, (zq, Co); 
since the evolution of the system of interest (z) is coupled to that of the heat bath {Q. 

In Ref. [0, it was shown that, if one started with a distribution of initial conditions in 
the extended phase space given by 

/(z, C, 0) cx exp -[/?^fo(z) + C'rV2], (36) 

and if one then propagated an ensemble of trajectories from these initial conditions, under 
the Nose- Hoover (NH) equations of motion, as A was switched from to 1 over a time t^, 
and computed the work W for each trajectory, then Eq.|^ would hold identically for this 
ensemble. This was shown by inspection. 
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Our purpose here is to use the central result of the present paper to establish simple 
criteria for determining whether or not Eq.|^ is valid for a particular implementation of IMD. 
In general, the heat bath is represented by N variables, Ci to (n. Let y = (z, (^i, ■ ■ ■ , (n) 
denote a point in the {2D + A^)-dimensional extended phase space, where D is the number 
of degrees of freedom of the system of interest. Then the thermostatting scheme in question 
is defined by a set of deterministic equations of motion in this space: 

y = KA(y), (37) 

of which Eq.^ is an example. An ensemble of trajectories evolving under these equations 
of motion may be described by a distribution /(y, t) which satisfies the continuity equation, 

|._iL.(K,/).«,/. (38) 

Now suppose that we can find a function q{Ci, ■ ■ ■ , (n) such that the distribution 

/f cx exp -P[H,iz) + g(Ci, ■ ■ ■ , Cn)] (39) 

is stationary under the evolution defined by Eq.|37|, when A is held fixed. That is, 

Rxfx = 0. (40) 

(For the NH equations, q = satisfies this condition.) We will allow q to depend on 

P, and any other constant parameters, but not on A. The distribution fxiy) may be viewed 
as the "canonical" distribution in the extended phase space, since it is invariant when A is 
held fixed. Now define a function 

nx{y)^H,{z) + q{Cu---XN), (41) 

which we may think of as an extended Hamiltonian. (This is not meant to imply that 
Tix generates Eq.^ as a "real" Hamiltonian would; in general those equations are non- 



Hamiltonian.) By Eq.iO 



^Aexp-/57^A(y) = 0. (42) 

Now forget for a moment the division between system of interest and heat bath, and treat 
the entire extended phase space as a phase space for some system of interest, for which a 
parameter-dependent energy function, Hxiy), is defined. The evolution given by Eq.|3^ then 
satisfies the two conditions listed in Section ||: (1) it is Markovian (Eq. |38|) , and (2) it satisfies 
detailed balance (Eq.^2|). Thus, our central result, exp—^W = exp — /?AF, is identically 



true for an ensemble of trajectories evolving under Eq.^ from an initial distribution /q oc 
exp —(3T-Cq, provided we replace H\{z) by l-t\{y) in computing W and AF. However, it is 
easily verified that we will get exactly the same values for both W and AF using H\{y), 
as we would obtain with Hx{z). In the case of W, this is because only the first term of H\ 
(namely, Hx) depends on A; for AF, it follows from the fact that the partition function for 
the extended Hamiltonian factorizes into an integral over the z variables and an integral 
over the ( variables, and only the former depends on A. (See the definitions of work and 
free energy, Eqs.^ and ^) These considerations lead to the following simple conclusion. 
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Suppose we have a system described by a Hamiltonian Hx{z), and we wish to compute 
the free energy difference AF = Fi — Fq. Suppose furthermore that we are given a scheme 
for IMD; that is, we have a set of "heat bath" variables (Ci, ■ ' ' ; Ca^); along with equations 
of motion in the extended phase space, as per Eq.^. Then, if a function q of the heat bath 
variables can be found such that the distribution oc exp —(3{H\ + q) is stationary under 
the equations of motion (with A fixed), then we can compute AF as follows. Let an ensemble 
of trajectories evolve from an initial distribution /o oc exp —(3{Hq + g), as A is switched from 
to 1 over a finite switching time tg- Compute the work W for each trajectory, as per Eq.^, 
and then compute the ensemble average of exp —[3W. This average will equal exp — /3AF. 



D. Monte Carlo evolution 

As a final example — again, motivated by computer simulations — let us consider the 
evolution of a Monte Carlo (MC) trajectory, as A is switched from to 1. In this case both 
the trajectory and the parameter A evolved in discrete steps, rather than continuously: 

z(t) Zo,Zi, ■ ■ ■ ,ZAr (43) 

A(t) ^ Ao,Ai,---,A;v ; K = n/N. (44) 

The initial point in phase space, zq, is sampled from a canonical distribution, with the 
value of A fixed at Aq = 0. One then imagines that A changes abruptly from Aq = to 
Ai = as a result, a quantity of work 5Wi = H\^{7.q) — Hxf^{zo) is performed on the 

system. Then, the system jumps to the next point in phase space, zi, generated from zq by 
a MC algorithm appropriate to the Hamiltonian Hx^ (see Appendix B). One then continues 
to alternate between discrete changes in A and discrete MC jumps in phase space, until 
the entire "trajectory" (zo,---,ZAr) is obtained, and the value of A is 1. The total work 
performed during this discrete switching process is: 

N N 

W = J2SWn = Y. [Hx„ (z„_i) - Hx„_, (z„_i)] . (45) 

n=l 71=1 

Note that "time" does not enter into this scheme. The quasi-static limit is obtained by 
letting the number of steps become arbitrarily large: N — >■ oo. 

Let us now imagine that we generate an infinite ensemble of MC trajectories, each of 
length A^. We do this by implementing the above procedure repeatedly, each time feeding in 
a different string of random numbers to generate the initial conditions zq, and the subsequent 
MC steps. Given such an ensemble, we compute the work W performed on each trajectory 
(Eq.^5]), and then the ensemble average of exp —(5W. As shown in Appendix B, this average 
will equal exp —(3AF. This should come as no surprise: a trajectory generated by the MC 
algorithm is a Markov chain, with detailed balance built into the individual steps. This 
evolution thus satisfies, in discretized form, the two assumptions required for the validity of 
Eq.|. 

Let us now consider the two limiting cases = 1 and A^ — * oo. The latter is the 
MC equivalent of the quasi-static limit. In this case our ensemble proceeds through a 
discrete, infinitesimally spaced sequence of canonical equilibrium distributions, and the work 
performed on each trajectory is AF, as per Eq.ffl. The result exp —j3W = exp —jSAF then 
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follows automatically. In the opposite limit, = 1, the work W performed on a particular 
"trajectory" is given by 



W = Hi{zo) - Ho{zo) = AH{zo) {N = 1). 



Eq.|3 then reduces to: 



(exp -pAH)o = exp -/?AF, 



(46) 



(47) 



where {■ ■ ■)o denotes a canonical average with respect to A = 0. This result is a well-known 
identity for the free energy difference AF; see Eq.ffS] below. 



IV. FREE ENERGY COMPUTATIONS 

While Eq.|]is interesting in its own right, it may additionally prove useful in the numerical 
computation of free energy differences. The field of free energy computations is decades old, 
with diverse applications, and a very large body of literature exists on the subject |T^. In 



this section, without attempting a survey of the field, we discuss a few points relevant to 
the possible application of Eq.|^ to free energy computations. These comments expand on 
ones made in Ref. |@. 

Most methods of computing free energy differences are variants of either the thermo- 
dynamic integration (TI) or thermodynamic perturbation (TP) methods. (For an exception 
to this statement, see the work of Holian, Posch, and Hoover [Q, where two new expres- 
sions for AF, based on time-integrated heat transfer, are derived within the framework of 
isothermal molecular dynamics.) TP is based on the identity |T3 



AF = -(3'Hn{exp-/3AH)o , {AH = H^- Hq), (48) 

where (■ ■ ■)\ denotes a canonical average with respect to a fixed value of A. Using a method 
such as Monte Carlo, one samples points in phase space from the canonical distribution 
corresponding to A = 0, and then one takes the average of exp —(3 AH over these points. 
In principle, the method is exact for N — > oo\ in practice, unless the canonical distributions 
corresponding to Hq and Hi overlap to a significant degree, the average of exp —(3 AH will be 
dominated by points in phase space which are visited extremely rarely during the canonical 
sampling, so numerical convergence with will be prohibitively slow. One way to get 
around this problem is to break up the A-interval [0, 1] into small sub-intervals, then use 
Eq.HB| to compute the free energy difference corresponding to each sub-interval. Other, more 



sophisticated methods of extracting the best efficiency from Eq.^ have been developed over 
the years |T^. 

TI is based on the identity 

The integral on the right may be evaluated by sampling Us points in phase space from each 
of M different canonical distributions (corresponding to equally spaced values of A from to 
1), for a grand total of A^ = UsM points sampled. The average of dH\/d\ is then computed 
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at each value of A, and from these M averages the integral is obtained. With the exception 
of a few very simple systems (e.g. ideal gases, harmonic oscillators), the standard way to 
obtain a canonical distribution is to first allow the system to "age" — that is, to relax to 
an equilibrium statistical state — under some MC or IMD scheme. In implementing the 
numerical evaluation of Jq dX {dHx/dX)x, it is often too time-consuming to age the system 
independently at each of the M selected values of A. Instead, the final point sampled at 
one A value may be used to generate the initial point sampled at the next value of A. Of 
course, this means that at each A (except A = 0) we are actually sampling from a slightly 
out-of-equilibrium distribution, which leads to a systematic error (in fact, an over-estimate) 
in the evaluation of the integral. 

A limiting case of this procedure arises when we take = 1, i.e. exactly one point is 
sampled at each value of A. This is the slow growth method [0. As stressed by Reinhardt 
and coworkers [3,P^,^ , if we view the chain of points (zo,---,Z7v) thus generated as a 
trajectory evolving in phase space (as A is switched from to 1), then the integral appearing 
on the right side of Eq.^ represents the work W performed on the system over the course 
of the switching process. This picture is a compelling one, as it attaches a very physical 
interpretation to the numerical evaluation of Eq.^ instead of computing an integral, we 
are simulating the evolution of a thermostatted system, with the idea that, in the limit 
of infinitely slow switching, the work performed on the system will equal the free energy 
difference AF. For a finite switching rate, the above-mentioned systematic error inherent in 
the slow growth method is simply a manifestation of the inequality W > AF; see Eq.|] and 
Refs. p|, |20|pl| . This interpretation of free energy computations is referred to as adiahatic 
switching (or finite-time variation) , and indeed may be viewed as a separate method, distinct 
from thermodynamic integration. 

In the context of adiabatic switching, Eq.^says that, if we run an ensemble of finite-time 
(or finite- A^, for Monte Carlo) simulations of the switching process, using for instance one 
of the methods described in Section then the average of exp —j3W over this ensemble 
of simulations will equal exp —(3AF. Assuming perfect numerical accuracy and an infinite 
ensemble of simulations, this is an exact statement. Another way of putting it is as follows: 
for a single finite-time switching simulation, the value of exp —f3W provides an unbiased 
estimate of exp —(3AF. By contrast, the value of W gives a biased estimate of AF (i.e. 
W > AF). Indeed, the former statement implies the latter, in the same way that Eq.|48| 
implies the Gibbs-Bogoliubov-Feynman inequality [^: {AH)q > AF. Let us consider this 
for a moment. 

Given a real function of a real variable, y{x), it is easy to show that, if d'^y/dx^ > for 
all X, then 

T N 1 AT 

i=i 1=1 

for any set of points {a;i, ■ ■ ■ ,XAr}. If these points are the result of random sampling from 



some ensemble, then Eq.50, in the limit N oo, may be rewritten as 



y{x) > y{x), (51) 

where the overbar denotes an ensemble average. Applying this result to y{x) = expx, with 
X = —PW, we get 
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exp -(3W > exp -(3W. 



(52) 



This is identically true for any distribution p{W). We may now combine this result with 
Eq.| to get W > AF. 

Now, what if we perform a finite number, Ns, of identical switching simulations? Let Wi 
denote the work performed on the system during the i'th simulation, and let 

^^■5 1 = 1 

be the average over these values. We may view the WiS as numbers sampled randomly from 
a distribution p{W) satisfying Eq.^. Then the expectation value of W"" provides a rigorous 
upper bound on AF: 

{{W)) = W = J dWp{W) W > AF. (54) 

The double angular brackets, denoting expectation value, specifically mean an average over 
all possible sets of Ns simulations. Now, Eq.H suggests that, rather than W"', we consider 
the following quantity as our best estimate of AF: 

1 Ns 

= In — E -m- (55) 

* i=l 

For Ns = 1, and W"' are identical, and the expectation value of either is W. For 
Ns — > oo, by contrast, converges to AF, whereas converges to W. For intermediate 
values of Ns, the following inequality chain holds: 

AF < {{W)) < {{W)). (56) 

(Both inequalities are derived by combining Eqs.^ and ^ with the definitions of and 
W^.) In other words, as an estimate of AF, the "exponential average", W^, is statistically 
less biased than the ordinary average, W^, for A^s > 1- 

On the face of it, the last statement seems to imply that, if we perform more than one 
switching simulation, then we are better off using rather than as our best guess (or 
upper bound) for AF. In practice, however, Eq.^ may be subject to the same disease as the 
TP identity, Eq.^. Namely, if the values of W obtained from repetitions of the switching 
simulation typically differ from one another by much more than j3~^, then the average of 
exp —f3W will be dominated by values of W which are very rarely sampled Thus, the 
convergence of to AF, in the limit Ns oc, may be much slower than the convergence 
of to W , in the same limit. In other words, for a finite number of switching simulations, 
may be subject to considerably larger statistical fluctuations than W"", even though its 
systematic error (expectation value minus AF) is, by Eq.^, smaller. 

The preceding comments point to the following tentative conclusion. If one runs a set of 
Ns switching simulations, with the goal of computing AF, and if the spread in the values 
of work W obtained is not much larger than (3~^, then the exponential average defined 
by Eq.^ should provide a better estimate of (or tighter upper bound on) AF than the 
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ordinary average W^. This conclusion is supported by a calculation by John E. Hunter III, 
as described in Ref. [Q; see also the numerical illustrations in the following section. 

Of course, a more detailed study of the possible utility of Eq.|^ to free energy computations 
should be made. In particular, it is not ruled out that there exist methods around the 
limitation mentioned in the previous paragraphs . 



V. NUMERICAL RESULTS 

In this section we illustrate our central result with numerical experiments. The first four 
sets of simulations involve a harmonic oscillator whose natural frequency is switched from 
LJo = 1.0 to uji = 2.0. The evolution is implemented using, in turn, each of the four examples 
discussed in Section |T|. Then, we present results involving a more complicated system: a 
gas of interacting particles inside an externally pumped piston. All of these cases satisfy the 
condition discussed at the end of Section |IV|, namely, the spread in the values of W is not 
much greater than (Otherwise, the convergence of to AF, in the limit of many 

simulations, would be poor.) 

For the harmonic oscillator simulations, we take the Hamiltonian given by Eq.^ with 
ujx = 1.0 + A, and A(t) = t/tg as throughout this paper. Also, we take (3"^ = 1.5. Thus, the 
free energy difference is AF = ln(a;i/a;o) = 1.0397. 

In the first set of simulations, the oscillator is isolated (it evolves under Hamilton's equa- 
tions, with A time-dependent), although the initial conditions are sampled from a canonical 
ensemble corresponding to A = 0: 

/(a;,p,0) = ^exp[-/3(p2 + ^oV)/2]. (57) 

Five different values of the switching time were chosen: ts = 1.0, 3.0, 10.0, 30.0, and 100.0, 
and for each ts a total of Ng = 10^ simulations were carried out. Fig.^ shows the average 
value of work obtained at each switching time, W"", as well as the "exponential average", 
W^. (See Eqs.^ and ^.) Since Hamiltonian evolution satisfies detailed balance, but not 
thermalization (as defined in Section |ID, we do not expect the work performed in the limit 
ts ^ oo to equal the free energy difference AF. Rather, we expect W^o = [(i^i/i^o) ~ l]-E'o 
for a single trajectory (see Section pi A|) , and therefore, for a canonical distribution of initial 
energies. 



lim W = [(uji/ujo) - l]/9"^ = 1.5. (58) 



The values of W shown in Fig.^ are consistent with this expectation. Eq.^, meanwhile, 
predicts 



lim W"^ = AF (59) 



for any value of ts. Again, the numerical results are consistent with the prediction: the 
values of shown in Fig.|^ all fall very close to AF. 

For the second set of simulations, we added a frictional and a stochastic force, as described 



by Eqs.p9^ to 31, with Dp = 0.6 and /? =1.5. The evolution now represents that of an 
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oscillator coupled to a heat bath. The stochastic force was implemented by generating a 
random momentum kick (sampled from a Gaussian distribution) at each time step, dt = 0.01, 
in the numerical integration. As with the Hamiltonian evolution, 10^ simulations were 
performed at each of the five switching times, and the results for W"' and are plotted in 
Fig.^. Here we do expect that the work W will approach AF = 1.0397 as ts oo, and the 
results for W"" support this. At the same time, the exponential average falls very close 
to AF for each switching time, as predicted by Eq.|^. 

The next simulations again involved a thermostatted harmonic oscillator, only this time 
isothermal molecular dynamics was used to implement the coupling to the heat bath. The 
particular IMD scheme used was developed by Hoover and Holian [^], and involves two 
heat bath variables, ( and ^. The equations of motion in the extended phase space are 

X = p (60a) 

p = -iulx-Cp-P^P^ (60b) 

C = iPp' - 1)/t' (60c) 

e = iPY - W)lr\ (60d) 

where r is a relaxation time whose value was set to unity. A total of Ng = 10^ switching 
simulations were performed, with a swtiching time tg = 1.0. At the start of each simulation, 
initial conditions were sampled from the distribution 

f{x,p, C, e, 0) cx exp[-/?(p2 + ujy)/2 - r'iC' + e')/2]. (61) 

It is easily verified by inspection that, if the value of A were held fixed at 0, then this 
distribution would be invariant under Eq.^. This is therefore the "canonical distribution" 
corresponding to this IMD scheme, and the function q defined in Section [HI C| is given by: 



qiC,0 = r\e + e)/2p. (62) 

This set of simulations was used to illustrate Eq.|2y, evaluated at t = tg. (See the 
discussion at the end of Section ||.) Fig.§ is a scatter plot showing the distribution of final 
values in phase space, {x{ts),pits)), for the 10^ trajectories. Fig.|| shows several contour 
lines of this distribution, after smearing each point with a Gaussian of variance 0.04 in both 
the X and p directions. Thus, the lines shown are actually contours of the function 

1 

f{x,p) = —Y,6,[x-x,{Q]5,{p-pM] , e = 0.04, (63) 

^^■5 i = l 

where is a normalized Gaussian of variance e, and {xi{t),pi{t)) gives the phase space 
evolution of the z'th trajectory in the ensemble of simulations. As can be seen from both 
figures, this distribution does not correspond to the canonical distribution for Hi. Indeed, 
its skewness illustrates the "lag" which develops between the evolving phase space density 
and the instantaneous canonical distribution. 

Next, the solid lines in Fig.^ show contours of the function g{x,p,ts) defined in Section 
^ obtained from the same set of simulations. Again, Gaussian smoothing was used, so the 
solid lines are contours of the function 
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g{x,p,ts) = -^^its)]^P-Pi{ts)]exp-(3Wi , e = 0.04, (64) 

where Wi is the total work performed on the system during the i'th simulation. The dashed 
lines in Fig.|| show the corresponding contours of the predicted "mass density" g{x,p,ts) 
(Eq.^Op, with the same Gaussian smoothing function folded in. The agreement between 
the two sets of contours is very good. This shows that, indeed, when one assigns a weight 
exp —(3Wi to each of the points in the scatter plot, Fig.^ then the resulting weighted dis- 
tribution is canonical, in the sense of Eq.|20|. 

In the final set of simulations involving the harmonic oscillator, we used Monte Carlo 
evolution, with the Metropolis algorithm. Here, the duration of a simulation is characterized 
by the number of MC steps, N , rather than by a switching time t^. Ten different values of 

were considered, N = 5, 10, 20, 50, ■ ■ ■ 5000, and for each a total of 10^ simulations were 
performed. Fig.^ shows W"' and for each value of A^; as before, the results agree nicely 
with Eqs.|T] to ^. Fig.|^ shows Pn(W) — the distribution of values of W obtained from the 
10^ simulations — for each of the ten values of A^. Although the distributions pN are quite 
different, the integral / dW Pn{W) exp —[3W (i.e. exp —[3W^) is independent of A^, as shown 
by the values of in Fig.|^. 

As a final example, we take a system more complicated than a harmonic oscillator, 
namely, a gas of rip = 50 interacting particles inside a piston which is taken through one 
cycle of pumping. Specifically, the particles are confined within a two-dimensional box with 
hard walls, whose initial dimensions are 1x1; over the course of the switching process, one 
of the wall first moves inward until the area enclosed by the box is three-quarters of its 
initial value, and then back out again. This pumping of the piston is cosinusoidal: the x— 
and y— dimensions of the box are given by 

^.(A) = 1.0 (65) 
Ly{\) = 0.875 + 0.125 cos(27rA), (66) 

where X = t/tg and the total switching time is = 10.0. 

In their interactions with one another and with the walls, the particles act as hard disks 
of radius R = .005; between collisions each particle moves freely. (Thus, over a switching 
time ts = 10.0, a typical particle suffers several collisions with other particles.) Work is 
performed on the gas each time a particle bounces off the moving wall. 

Molecular dynamics was used for the evolution, i.e. continuous trajectories for the par- 
ticles were computed as functions of time. However, at each time step in the integration of 
the equations of motion, a single particle was randomly selected, and a random "kick" — 
a discrete change in the momentum of the chosen particle — was generated. The kick was 
then either accepted or rejected according to the Metropolis algorithm |2^, corresponding 



to a temperature = 0.5. This simulates coupling to a heat reservoir: if all the walls of 
the box were fixed, the gas would relax to thermal equilibrium from any initial conditions. 

Fig.^ shows results obtained from Ng = 10^ such simulations. For each simulation, the 
initial conditions of the gas were chosen from a canonical ensemble (/3~^ = 0.5), and the 
work performed on the g function of time, was computed. The horizontal axis shows 

time. The solid line gives the work performed on the gas up to time t — i.e. the "work 
accumulated", w{t) — averaged over all A^^ simulations. Let w"'{t) = (1/A^s) ^^^(t) 
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denote this average, where Wi{t) is the work accumulated during the z'th simulation. The 
dashed line gives the "exponential average" 



1 Ns 

w%t) ^ -p-' In — E exp -Pw,{t), (67) 

of the work accumulated. Both w"'{t) and w^{t) were computed from the same set of 10^ 
trajectories. 

Since the piston returns to its initial position at the end of the switching process, the 
final free energy difference is zero: AF = Fi — Fq = 0. We see that the dashed line indeed 
returns to zero ai t = tg, with very good accuracy. By contrast, the average work performed 
on the system (upper line) ends at W"' = w"'{ts) = 1.534. This represents dissipated energy: 
the gas "heats up" when pumped at a finite rate. 

At intermediate times, we expect 

w^{t)=F^-Fo , A = A(t) (68) 



(in the limit of infinitely many switching simulations), by Eq.|2^. If the gas were truly ideal, 
then the free energy difference would be given by: 

Fx-Fo = Up (3-^ HAo/Ax) (ideal gas) , (69) 

where A\ denotes the area enclosed by the box. However, since the particles do interact 
with one another, as hard disks, this expression for Fx — Fq is not exact. Nevertheless, the 
size of each particle is small enough {R = .005) that Eq.^ ought to represent an excellent 
approximation. The dotted line in Fig.^ shows this approximation to Fx — Fq. This line 
is very close to the exponential average (dashed line), in confirmation of our central result, 
in the form given by Eq.|22|. Note that the dotted line represents the work that would have 
been performed on the gas, in the limit of infinitely slow switching (t^ — > oo). Thus, in the 
dahsed line, we have effectively extracted this quasi-static behavior, from an ensemble of 
finite-time switching simulations! 



VI. SUMMARY AND DISCUSSION 

The central goal of this paper has been to establish the (exact!) validity of the result 
exp —(5W = exp —f3AF, within the framework of the master equation approach. This 
result is unusual, in that it expresses in the form of an equality (rather than an inequality, 
e.g. Eq.^ the relationship between the work W performed on an out-of-equilibrium system 
(more precisely, on an ensemble of systems driven out of equilibrium by varying an external 
parameter), and the free energy difference AF between two equilibrium states of the system. 
A few comments are now in order. 

In classical statistical mechanics, the equilibrium "state" of a system is described by a 
canonical distribution in phase space: /*"(z) = exp —j3H{z). Its free energy is then 
given by 

F = (H) - = -(3~HnZ, (70) 
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where (H) = J H and S = —ks J In/*^. Thus the free energy F (hke the entropy S) is 
a quantity associated with a statistical ensemble of microscopic states of the system. When 
the system depends on some external parameter, then so does the canonical ensemble, as 
in turn does the free energy. The quantity AF of central interest throughout this paper 
has been the free energy difference (at constant temperature) between two such equilibrium 
ensembles, /f=o(z) and /f=i(z). 

In deriving Eq.^ we introduced a time-dependent phase space density /(z,t), describing 
the evolution of our ensemble of trajectories. Note that, although the initial phase space 
density coincides with the canonical distribution used to compute the free energy Fq, 

/(z,0) = /to(z), (71) 

the final density does not (typically) coincide with the distribution for which Fi is defined: 

/(z,t,)^/t,(z). (72) 

This is due to the lag which develops between the ensemble of trajectories and an instan- 
taneous canonical distribution in phase space. Thus, AF is not the free energy difference 



between the initial and final states of the system ^Tj , but rather — as stated in the previous 
paragraph — between two different canonical ensembles, /£^q and fx=i, only the former of 
which reflects the actual distribution of microscopic states of the system at any time during 
the switching process. 

Another way of putting this is as follows. Suppose we are interested in the free energy 
difference between two equilibrium statistical states of a system, A and B (corresponding 
to and fx=i, respectively). Ordinarily, we would compute or measure AF by reversibly 
carrying the system from A to B, i.e. by switching A infinitely slowly. Eq.^ tells us that, 
even if we switch irreversibly, so that the system ends up in some nonequilibrium statistical 
state B*, we can still extract AF{= F^ — F^) from an ensemble of such measurements. 

[Of course, if we are dealing with a system which satisfies the thermalization assumption, 
then, at the end of the switching process, we can always, at no cost in work, hold the value of 
A fixed and allow our ensemble to relax, for an additional time tj-ei, to thermal equilibrium: 
f{z,ts + trei) = /a=i(z). In this case the AF which appears on the right side of Eq.| does 
equal the free energy difference between the initial and final statistical states of the system.] 

Eq.^ was derived, as an identity, under the assumptions of Markovian evolution and 
detailed balance, as spelled out in Section |. This derivation is complementary to the one 
presented in Ref. [Q, in which the degrees of freedom of the heat reservoir were treated 
explicitly. Neither of the assumptions of Section | was assumed in Ref. 0] , but the coupling 
between system and reservoir was taken to be weak, so the result there was an approximate 
one, with small corrections expected from the small but finite interaction Hamiltonian. 
Of course, in a real physical system, neither the Markov assumption nor detailed balance 
will be met exactly, so the derivation presented herein is strictly speaking valid only for a 
particular class of models of physical reality. Nevertheless, because the result is exact for 
these models, and because the Markov and detailed balance assumptions are often very good 
approximations for physical systems, the result is a useful one. Furthermore, as illustrated 
in Section models of thermostatted systems which are commonly used in theoretical 
and numerical studies do satisfy these assumptions; Eq.^ is therefore exactly valid for these 
models. 
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It would be very interesting, of course, to find a pliysical system on which a laboratory 
(as opposed to numerical!) experiment testing the validity of Eq.^] would be feasible. As 
mentioned in Ref. such a system would almost certainly have to be micro-, or at most 
meso-, scopic in size. 

Finally, from both a theoretical and a computational point of view, it would be worth- 
while to consider possible extensions or generalizations of Eq.|[ In particular, are analogous 
results valid for ensembles other than the canonical ensemble (fixed N, V, T) considered 
here, e.g. microcanonical, grand canonical, isothermal-isobaric, etc? Presumably, the role of 
the Helmholtz free energy F would then be played by other thermodynamic potentials, for 
instance the Gibbs free energy in the case of the isothermal-isobaric ensemble. 
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APPENDIX A 

In this Appendix we present a derivation of Eq.^ different from the one given in Section 

For a given stochastic trajectory z(t), the work W is given by a path integral along that 
trajectory (Eq.^. We are interested in evaluating the average of exp —jSW over an ensemble 
of trajectories, obtained by sampling initial conditions from a canonical ensemble and then 
evolving stochastically from each of these initial conditions. The quantity exp —(3W thus 
constitutes a "sum over all paths", with each path z(t) in our ensemble weighted by the 
factor exp —(3W. We may write this as 

exp -^W = J dmV[z{t)] exp -(3W, (73) 

where dm denotes a measure in the space of paths z(t); V[z{t)] denotes the probability 
density (with respect to this measure) of choosing z(t) by sampling randomly from the 
ensemble of trajectories; and W = W[z(t)] (as per Eq.|^). 

Let us now divide the time interval [0,ts] into N time steps of duration 6t = tg/N, and 
let us denote a particular trajectory z{t) by its phase space locations z„ = z(t„) at times 
tn = n6t,0<n<N. Thus, 

Z(t) (zo,Zi, ■ ■ ■ ,ZAr). (74) 

The limit N ^ oo (with tg fixed) is implied. Choosing a Euclidean measure in path space, 

dm = J dzo J dzi ■ ■ ■ j dz^, (75) 
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the probability density for a particular path is 



V[z{t)] = p(zo)P,^:(zo|zi)Pa1(zi|z2) ■ ■ ■P,^^(z^-i|z^.). 



5t, 



(76) 



Here, p(zo) = Zq^ exp —PHq(zq) is the probability distribution for the initial condition zq; 
P^*(z'|z) is the transition probability from z' to z (in time St) as a function of A; and 
Xn = n/N. It is the Markov assumption which allows this factorization. The work W may 
be expressed as 



N 



(77) 



n=l 



where 5Hn = Hx^ — H\^^. [In writing Eqs.|7g and |73 we implicitly assume that \{t) evolves 
in discrete steps 5\ = 1/N which occur at times to, ^i, ■ ' ' , ^iv-i- This "staircase" evolution 
becomes X{t) = t/tg in the limit — » oo.] 
Combining Eqs.l73| to 17^ we arrive at 



N 



exp —PW = 

Ln=0 ■ 

Let us now introduce 

rM-i 



p(zo)e-^^^^(^°)p,^;(zo|zi) ■ ■ ■e-^^^-(^--)p,^^(z^._i|z 



N) 



(78) 



9m{z) 



n 



n=0 



1Z„ 



p(zo)e-'^^^^(^")pi*(zo|zO ■ ■ ■e-^^^-(^«-)p,^^(zM-i|z), 



(79) 



where 1 < M < A^. This is the discretized version of the function g{z,t) introduced in the 
main body of the text: 



5'Af(z) = giz,t 



M) 



(80) 



In particular, note that exp —(3W = J dzgi^{z). This set of functions qm satisfy the recursion 
relation 



(?M+l(z) = / t/ZM^M(zM)e-'^'^-^+^(^")P,t^^(zM|z). 

Now, to first order in 5t, we have 

e-^'''''+'^^''^ = 1- P5Hm+iM 
P\m+i(^m, z) = 6{zm - z) + 5tPAM+i(zA/, z). 

Combining this with our recursion relation gives (to leading order) 

^[gAf+i{z) -gM{z)] = -pgM{z) 5HM+i{z)/6t + j dzAf 5-^^ (za/) Pam+i(za/, z), 
which becomes Eq.|l6| in the limit oo. 



(81) 



(82) 
(83) 



(84) 
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APPENDIX B 



Here we prove the assertion (made in Section pn| ) that Eq.^ is identically true when the 
switching process is carried out using the Monte Carlo method. Some of the steps in the 
proof will be similar to those in Appendix A, but the assumption N oo will not be made 
here. 

As mentioned in Section |T|, a trajectory (zq, • • • , z^v) is obtained by alternating discrete 
changes in the value of A with random jumps in phase space generated by the MC algorithm. 
This algorithm — parametrized by the value of A — takes as input a point z, and outputs 
a point z'. Let Px{z\z') denote the probability of generating an output z' from an input z, 
for a given value of A. Detailed balance is built into the algorithm: 

J rfze-^^^(^) Px{z\z') = e-^^^(^'), (85) 

for any A. (This may be accomplished by, e.g., the Metropolis method |^6|.) Thus, a 
canonical distribution of inputs z gives a canonical distribution of outputs z'. 

The probability of obtaining a particular trajectory (zq, ■ ■ ■ , z/yr) over the course of the 
entire switching process is then 

P(zo, ■ ■ - z^v) = ^e-^^«(^«) Pa,(zo|zi) • ■ ■ PxA^n-i\zn)- (86) 
Combining this with Eq.^ for the work, we get 

exp-(3W= n — e-^^°(^»)e-''^^^(^°)PAi(zo|zi)---e-^^^-(^--)PA^(z;v-i|z^), (87) 

.n=0-' J ^0 

where 5Hn = Hx„ — Hx„_i- Now notice that exp — /9ifo(zo) can be combined with 
exp —P6Hi{zq) to give exp — /3ifAi(zo)- The only other factor in the integrand which de- 
pends on Zq is Pxi(zo|zi). Performing the integral / dzQ, we get 

j o?zoP(zo|zi) exp -/3ifAi(zo) = exp -pHx^izi), (88) 

using Eq.^ This takes care of the first of the + 1 integrals appearing in Eqj8^. We now 
repeat this process, first combining exp —(3Hx^{zi) (obtained from the dz^ integration) with 
exp — /35if2(zi) to get exp — /3ifA2(zi), then integrating over zi, and so forth. At the end of 
this process of "rolling up" the factors and integrating, we are left with 

/■ 1 Z 

exp -(3W = dzN — exp -(3Hi{zn) = ^ = exp -(3AF. (89) 

J Zq Zq 

Q.E.D. 

Note that Eq.^ for g{z,t), derived within the framework of continuous-time evolution, 
also has a Monte Carlo counterpart. Namely, for 1 < M < iV, let us define 

M-l „ 

9m{z) = n J d'^^ 'Pm{zo, ■ ■ ■ , zm-1, z) exp -f3wM, (90) 

n=0 
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where Vm gives the probabihty of samphng a particular sequence of phase space points in 
the first M Monte Carlo steps, and wm is the work accumulated during those steps. Thus, 
in terms of our ensemble of MC trajectories, guiT) is the weighted phase space density 
after M steps, where the weight assigned to each trajectory is exp —[3wm- Then, writing an 
explicit expression for Vm in the form given by Eq.^ and rolling up factors and integrating 
as above, it follows easily that 




(91) 







■0 
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FIGURES 



FIG. 1. Distribution of values of work, p(W, tg), performed during an ensemble of independent 
switching measurements at a given switching time ts- The vertical line represents a delta function 
at W = AF, and corresponds to tg oo; in that limit, the work performed during a single 
switching process is exactly equal to AF. The smooth distribution represent p{W,ts) for a finite 
value of ts- In this case the ensemble average work exceeds the free energy difference, W > AF, 
since energy is dissipated in a finite-time (irreversible) process. 

FIG. 2. Simulations of an isolated harmonic oscillator whose natural frequency is switched 
from 0^0 = 1.0 to = 2.0 over a switching time tg- At each of five values of ts, 10^ simulations 
were carried out. The upper and lower sets of points show the ordinary averages (W"") and the 
exponential averages (W^) of the work, respectively. The dashed line is at VF = 1.5, the dotted 
line at IF = AF = 1.0397. 

FIG. 3. Same as Fig.^, except that the harmonic oscillator is now subject to a frictional and 
a stochastic force, as per Eq.^. The dashed line gives the free energy difference, AF = 1.0397. 

FIG. 4. In these simulations, the harmonic oscillator is thermostatted with the IMD scheme 
described in the text. 10^ simulations at tg = 1-0 were performed, and the dots in this figure show 
the final locations in phase space of these trajectories. 

FIG. 5. Contour plot of the distribution f{x,p, tg), constructed from the data shown in Fig.^, 
with Gaussian smoothing. 

FIG. 6. The solid line shows contours of the function g{x,p,ts), constructed from the data 
shown in Fig.^ the dashed line shows contours of the theoretical prediction for g{x,p,ts) (Eq.pO|). 
Both are smoothed with Gaussians. 

FIG. 7. Similar to FigJ^, except here the evolution of the thermostatted oscillator is imple- 
mented using Monte Carlo, rather than Langevin evolution. The duration of a simulation is now 
characterized by the number of MC steps, N, rather than a switching time tg. For each of ten 
values of N, 10^ simulations were carried out. 

FIG. 8. For each of the ten sets of MC simulations (see FigJ^, the distribution of values of 
work, /Oat (IF) was obtained. This figure shows these ten distributions, from N = 5 (lowest peak) 
to = 5000 (highest). (Although the peak moves toward the right with increasing A^, the actual 
average work performed goes down; see the values of in Fig.^.) 



28 



FIG. 9. Simulations of a gas of interacting hard disks inside a piston which goes through 
one cycle of pumping (first in, then out), over a switching time tg- The evolution is MD, but MC 
"kicks" are included to provide a thermostat. A total of 10^ simulations were performed. The solid 
line gives the average, and the dashed line the exponential average, of the work as a function of 
time, w"'{t) and w^{t), respectively. The dotted line gives the theoretical prediction for Fx — Fq, 
Eq.|69|, with A = X{t), for the case of an ideal gas. 
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